# Sesión 2: Resolución de ecuaciones no lineales

kill(all);

## El método de bisección.

Queremos resolver la ecuación 2*x+1+sin(x)=0

f(x):=2*x+1+sin(x)$

Utilizamos una gráfica para determinar un intervalo adecuado

plot2d([f(x)], [x,-1,0],
 [plot_format, gnuplot],
 [gnuplot_preamble, "set zeroaxis;"])$

Primer intento:

a:-1.0$b:0.0$

for k:1 thru 20 do 
(
  c:(a+b)/2,
  if f(a)*f(c)<0 then b:c else a:c 
);

c;

Ahora mostramos la sucesión de intervalos generados:

a:-1.0$
b:0.0$
for k:1 thru 20 do 
(
  c:(a+b)/2,
  if f(a)*f(c)<0 then b:c else a:c,
  print(k,a,b,c)
)$

Mostramos las sucesivas aproximaciones junto con los valores de f en ellas:

a:-1.0$
b:0.0$
for k:1 thru 20 do 
(
  c:(a+b)/2,
  if f(a)*f(c)<0 then b:c else a:c,
  print(c,f(c))
)$

Disponemos la misma salida con formato:

a:-1.0$
b:0.0$
for k:1 thru 20 do 
(
  c:(a+b)/2,
  if f(a)*f(c)<0 then b:c else a:c,
  printf(true,"~9,16f ~t ~9,16f ~%",c,f(c))
)$

Introducimos criterios de parada: número máximo de iteraciones, precisión y tolerancia

a:-1.0$
b:0.0$
max_iter:50$
tol:1.0*10^(-5)$
prec:1.0*10^(-3)$
for k:1 thru max_iter do 
(
  c:(a+b)/2,
  if abs(b-a)<tol then return(c),
  if abs(f(c))<prec then return(c),
  if f(a)*f(c)<0 then b:c else a:c,
  printf(true,"~9,16f ~t ~9,16f ~%",c,f(c))
)$

## Una rutina para el método de bisección.

Versión elemental:

biseccion(f, a1, b1) := block( [a, b, c],
    a : float(a1),
    b : float(b1),
    for n:1 thru 15 step 1 do (
        c: (a+b)/2,
        if f(a)*f(c) < 0 
            then 
                b : c 
            else 
                a : c
        ),
    c
    )$

biseccion(f(x),-1,0);

biseccion(f,-1,0);

Controlamos la viabilidad del intervalo inicial:

biseccion(f, a1, b1) := block( [a, b, c],
    a : float(a1),
    b : float(b1),
    if f(a)*f(b) > 0 
        then 
            return ("Intervalo no adecuado"), 
    for n:1 thru 15 step 1 do (
        c: (a+b)/2,
        if f(a)*f(c)<0 
            then
                b:c 
            else 
                a:c 
        ),
        c 
    ) $

biseccion(f,-1,0);

biseccion(f,0,1);

Controlamos que la raíz no se alcance ni en los extremos del intervalo inicial ni en ninguno de los puntos medios

biseccion(f,a1,b1):= block( [a, b, c, m, maxiter:300, tol:10^(-15), prec:10^-10],
    a : float(a1),
    b : float(b1),
    c : (a+b)/2,
    if abs(f(a)) < prec then return([0,{a}]),
    if abs(f(b)) < prec then return([0,{b}]),
    if f(a)*f(b) > 0 then return ("Intervalo no adecuado"),
    if abs(f(c)) < prec then return([1,{c}]),
    for n : 1 while ( n <= maxiter and abs(b-a) > tol and abs(f(c)) > prec) step 1 do (
        m : n,
        c : (a+b)/2,
        if f(a)*f(c) < 0 then b : c else a : c 
        ),
    [m,{c}] 
) $

biseccion(f,-1,0);

## Una función para el método de Newton.

Para programar el método necesitamos implementar la función derivada.

kill(all);

f(x):=x^3-5$

diff(f(x),x);

fd(x):=diff(f(x),x);

Podemos comprobar que de este forma no obtenemos una expresión evaluable: 

fd(3);

Hay varias formas de arreglar este problema:

fd2(x):= ''(diff(f(x),x));

fd2(3);

define(fd3(x),diff(f(x),x));

fd3(3);

También podemos utilizar el comando subst() para evaluar la derivada, que es lo que usamos en esta rutina 
que implementa el método de Newton:

newton(expr,var,punto_inicial):=block(
[x0:punto_inicial,x1,expd,iter,maximo_iter:15,precision:10^(-6),tolerancia:10^(-7)],
expd:diff(expr,var),
for k:1 thru maximo_iter do 
   (
    iter:k,
    if abs(subst(x0,var,expd))<precision then return("Error"),
    x1:x0-subst(x0,var,expr)/subst(x0,var,expd),
    print(k,x1,subst(x1,var,expr)),
    if abs(subst(x0,var,expr))<precision then return(x1),
    if abs(x1-x0)<tolerancia then return(x1),
    x0:x1
),
if iter=maximo_iter then
   print("Tras", maximo_iter,"iteraciones, el valor de la raíz es",x1),
x1
)$

newton(x^2-2,x,1.4);

Como alternativa podemos resolver el problema declarando funciones locales al módulo:

newton(expr,var,punto_inicial):=block(
[x,x0:punto_inicial,x1,expd,iter,maximo_iter:15,precision:10^(-6),tolerancia:10^(-7)],
local(f,fp),
define(f(x),subst(x,var,expr)),
define(fp(x),subst(x,var,diff(expr,var))),
for k:1 thru maximo_iter do 
   (
    iter:k,
    if abs(fp(x0))<precision then return("Error"),
    x1:x0-f(x0)/fp(x0),
    print(k,x1,f(x1)),
    if abs(f(x0))<precision then return(x1),
    if abs(x1-x0)<tolerancia then return(x1),
    x0:x1
),
if iter=maximo_iter then
   print("Tras", maximo_iter,"iteraciones, el valor de la raíz es",x1),
x1
)$

newton(x^2-2,x,1.4);